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ABSTRACT 

We examine the evolution of highly eccentric, planet-crossing orbits in the restricted 
three-body problem (Sun, planet, comet). We construct a simple Keplerian map in which 
the comet energy changes instantaneously at perihelion, by an amount depending only on 
the azimuthal angle between the planet and the comet at the time of perihelion passage. 

This approximate but very fast mapping allow us to explore the evolution of large ensembles 
of long-period comets. We compare our results on comet evolution with those given by the 
diffusion approximation and by direct integration of comet orbits. We hnd that at long 
times the number of surviving comets is determined by resonance sticking rather than a 
random walk. 


Subject headings: comets: dynamics - orbits: resonances - methods: numerical 


1. Introduction 

The orbits of small bodies in the outer solar system, such as planetesimals or comets, evolve 
through gravitational interactions with the giant planets. This process determines the structure of the 
Oort Cloud, the rate of depletion of the Kuiper belt, the total population and spatial distribution of the 
Centaurs, and the dynamical properties of observed comets. Understanding the nature of planet-induced 
orbital evolution is therefore important for a wide range of solar system problems. 

A particularly important issue is the evolution of highly eccentric orbits whose perihelia penetrate 
the region of the outer planets (e.g. long-period comets). In this case the evolution can be treated 
approximately as a one-dimensional random walk in energy; the energy change per perihelion passage is 
assumed to be a Gaussian random variable and the rms energy change is assumed to be small compared 
to the typical orbital energy. This “diffusion approximation” was introduced by Oort (1950) and has 
been employed by many authors, including Whipple (1962), Weissman (1978), Yabushita (1980), Bailey 
(1984), and Duncan et ah (1987). A closely related approach is to treat the evolution as a random 
walk; thus the energy change at perihelion is assumed to be a random variable but it is not assumed to 
be small compared to the orbital energy (Kendall 1961, Arnold 1965, Opik 1976, Everhart 1977, 1979, 
Yabushita 1979, Froeschle and Rickman 1980). The diffusion or random-walk approximation provides 


2 


considerable insight into the evolntion of comet orbits, bnt neglects potentially important effects snch 
as secnlar evolntion, resonances, etc. 

Direct nnmerical integrations provide the most accnrate way to analyze orbital evolntion. The 
evolntion of long-period comets nnder the gravitational inflnence of the giant planets has been examined 
by Wiegert and Tremaine (1998), and extensive integrations of short-period comets and Kniper-belt 
objects have been carried ont by Holman and Wisdom (1992) and by Dnncan and Levison (1997). 
Nnmerical integrations reqnire more work, which can be prohibitive, and may yield less insight than 
simpler methods. Thus it remains worthwhile to investigate approximate models for planet-induced 
orbital evolution. 

Clearly the random-walk approximation fails for orbits with semimajor axes comparable to those of 
the outer planets, since in this case the energy changes are correlated at successive perihelion passages. 
A simple approach that illuminates this transition is to replace the random energy change at perihelion 
in the random-walk approximation with an energy change that depends on the azimuthal phase of 
the planet at perihelion passage; thus the energy at successive perihelion passages is determined by a 
deterministic twist map of the form 


En+l — En + F^'ll^n), 

T y (1) 

where E is the specihc energy, Dp is the planet’s angular speed, P{E) oc \E\~^P is the comet’s orbital 
period, ip is the azimuthal angle between the planet and the comet perihelion at the time of perihelion 
passage, and F represents the energy changes caused by the planet at perihelion passage. This map was 
examined by Petrosky (1986) using a simple—but unrealistic—sinusoidal form for F] Petrosky called 
(|^) the Keplerian map, a term that we shall also adopt. Sagdeev and Zaslavsky (1987) derived the 
map independently, evaluating F in the simple but uninteresting limit that the perihelion distance was 
much larger than the planet’s semimajor axis (they also neglected the indirect term in the gravitational 
potential from the planet). Chirikov and Vecheslavov (1989) used the map to examine the long-term 
evolution of Halley’s comet. 

The aim of this paper is to provide a fairly rigorous derivation of the Keplerian map, to discuss its 
generalizations, to provide better approximations for the “kick function” F than have been available 
so far, and to compare the predictions of the map both to numerical integrations and the diffusion 
approximation, to clarify the strengths and weaknesses of each approach. One of our principal results is 
that the number of comets that survive for many orbits is much larger than predicted by the diffusion 
or random-walk approximation, a result that we ascribe to “resonance sticking”. 

Mostly we shall examine the motion of comets in a planetary system containing the Sun and a 
single planet of mass Mp on a circular orbit. Without loss of generality we may choose the period 
and the semimajor axis of the planet to be unity. With this choice of units the planet’s mean motion 
Dp = 27r, G{Mq -|- Mp) = dvr^, P{E) = /{\E\^P) is the comet’s orbital period, and the specihc 

energy of the planet Ep = —27r^. We also write rup = Mp/Mq. 
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1.1. The diffusion approximation 


In this approximation the evolntion of the comet orbit is treated as a one-dimensional random walk 
in energy, the mean-sqnare energy change per perihelion passage D is assnmed to be independent of 
energy and is assumed to be small compared to the typical comet energy E. Then the number of 
comets bound to the solar system at time t with energy in the range [E, E + dE\ is n{E, t)dE, where 


dn 1T-, 

'm ^ [p{E) 


( 2 ) 


A more convenient form is 


dn _ 2 n 

dt L 


{-Ef^% 


where Do is independent of the planetary mass, 

nipDo — 25 / 2^3 


D = 


25/%3 




( 3 ) 

( 4 ) 


and the average (■) is taken over the phase of the planetary orbit at a single perihelion passage. 

The solution to equation depends on the boundary conditions. The simplest assumption is 
that all comets with positive energy escape, so n = 0 for E > 0. In this case the Green’s function 
corresponding to the initial distribution n{E,t = 0) = 6{E — Eq), Eq < 0, is (Yabushita 1980) 


n{E, t) 


2 |Eo|V 2 


4 

mlDot 


Eo\d‘^+ \E\P^ 
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mlDot 



( 5 ) 


where I 2 is a modihed Bessel function. The total number of surviving comets after time t is given by 
the depletion function 


/■LXJ j- , 

N(t) = / n(E, t)dE = 72 , 4|B„|'''y (mjBoi) 


( 6 ) 


where 7 is an incomplete gamma function. At large times N{t) —>■ 8\Eo\ j(rripDot) . Note that equations 
d^) and are normalized so that the initial number of comets N{t = 0) = 1. 


A more accurate approximation is that comets are lost if either > 0 or < Emin] the absorbing 
boundary at Emin arises because comets in tightly bound orbits may collide with the Sun, evaporate, or 
evolve much more quickly under the influence of more massive planets. In this case it is straightforward 
to solve equation numerically to hnd the depletion function. 


2. Mapping method 

We examine the orbit of a test particle (the comet) moving in the combined gravitational held of 
the Sun and a single planet on a circular orbit (the restricted three-body problem). We restrict our 
attention to comets with small inclinations, since we are mostly interested in scattering of comets in 


















4 


the protoplanetary disk. Thus the comet’s phase-space position is described by its semimajor axis a, 
eccentricity e, argument of perihelion zu, and mean anomaly £. Alternatively we may use canonical 
elements. The simplest of these are the coordinate-momentum pairs (£, K = [GMqo]^^'^) and {w,L), 
where L = [GM 0 a(l — is the specihc angular momentum. The azimuth of the planet at time t 

may be written 0 p(f) = Qpt + (j)o = 27it -|- 0 o- 

The equations of motion are described by a Hamiltonian 

Q 4 

H{K,L,£,w,t) = + (7) 

For our purposes it is more convenient to use mean anomaly instead of time as the independent variable. 
Then the equations of motion are still Hamiltonian if we use {vu, L) and (f, —E) as canonical coordinates 
and —K{L, E, zu, t, i) as the Hamiltonian, where K is defined implicitly hj E = H{K, L, i, zu, t) and E 
is the energy. In other words 

dL _ dK dE dK dzu dK dt _ dK 
'di~dw' 'di~~~m' 'di~~~dL' 


We shall now focus on long-period comets. These spend most of their time at distances much greater 
than the planet’s semimajor axis, where they travel on near-Keplerian orbits around the Sun-planet 
barycenter. Changes in the comet’s energy and other orbital elements are localized near perihelion, 
where its interactions with the planet are strongest. Therefore it makes sense to assume that the 
encounter with the planet lasts for only a short time near perihelion (£ = 0). Furthermore since rup <C 1 
the perturbation from the planet is weak. Thus we may write the negative of the Hamiltonian as 

47r2 

K{L, E, zu, t, 1) = ^)<^ 27 r(^) + 0(mp), (9) 


where 52Tr{x) = ~ 2,Tzn) is the periodic Dirac delta-function and k = — 47 r^(—/ Hid£. 

Furthermore, k can only depend on the time through the planet azimuth (j)p{t) and can only depend on 
the azimuthal angles (pp{t) and zu through the combination = (j)p{t) — zu = 27 zt + (pQ — zu. Thus we 
may write 


K{L, E, zxj, t, t) = 




-2E)V2 


-|- mptz{L, E, 27it + 4>o — xzj) 62 n{£)- 


Equations (||) and (|l0f) imply that the Jacobi constant 


( 10 ) 


J = E — 2tzL 


( 11 ) 


is conserved, a property that we inherit from the restricted three-body problem. Note also that k is 
conserved as the trajectory crosses i = 0. 

To proceed further we shall adopt a simple functional form for k{L, E, ip) that (we hope) captures 
the most important physics of the interaction. A natural hrst approximation is that k is independent of 
E and L, so that k = K{'ip) (we discuss more general forms in §^). The assumption that k is independent 
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of energy is reasonable for long-period comets, which all have near-parabolic orbits in the vicinity of the 
planets. The assumption that k is independent of angular momentum is less natural since the interaction 
of a long-period comet with the planets depends strongly on its perihelion distance; however, the same 
assumption is made in the diffusion approximation (eq. and conservation of the Jacobi constant (eq. 
0 ) ensures that variations in L are small so long as the comets remain long-period (i.e. \E\ -C 1). 

Hamilton’s equations (§) can now be integrated from i = 27m — 0 to 27r{n -|- 1) — 0: 


zu = const, 

T'n+l Eji -|- 7Tlpf{lpn)) 

Tfl 

Ln+l f y 

V^n+l '^n T ^n)y 


where the “kick function” /(V’) = —2TidK,{ijj)/dxp is independent of planet mass, and the map we have 
derived is essentially the Keplerian map (2). 

Thus we have arrived at a symplectic mapping in two dimensions {'ijj,E), which depends on the 
planet mass and the kick function. The angular momentum L is determined by the energy and the Jacobi 
constant J (eq. 0). The derivation of the map was motivated by long-period comets, whose orbital 
period is much longer than the planet orbital period {P{E)/P{Ep) = (27r^/|i?|)^/^ ^ 1); however, it 
should also provide a fair representation of the behaviour of orbits with shorter periods, |i?|/(27r^) ~ 1. 


Mappings such as (|^ can have trajectories of two kinds, regular and stochastic. The energy of 
regular orbits will vary only over a limited range, whereas stochastic orbits can wander over the entire 
stochastic region, and in particular they can reach escape energy > 0), which corresponds to loss of 
a comet from the solar system. The comet is also lost if L < 0, which corresponds to collision with the 
Sun and occurs if < J. 


2.1. The kick function 

To determine the kick function /('0) we numerically integrated the orbits of a set of comets on 
initially parabolic orbits through a single perihelion passage. The comets shared a common perihelion 
distance q and were distributed uniformly random in argument of perihelion and longitude of ascending 
node. The inclinations were chosen at random from a narrow Gaussian distribution with zero mean and 
dispersion 0.1 radians. Then we calculated the energy change AE = Efinai — Einiuai in the barycentric 
frame and averaged AE over all comets with a given value of 'ipy which is the azimuthal angle between 
the comet and planet at the instant of perihelion passage. The kick function is the average energy 
change normalized by the planet mass, /(V’) = {AE)/mp. Figure 1(a) shows a scatter plot of the 
normalized energy change AE{'ip)/mp versus V’ for 30,000 passages with perihelion distance q = 0.5. 
The Figure also shows the corresponding kick function {AE)/mp, averaged over 2000 comets at each of 
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600 values of 'ijj. As we should expect, the kick function is approximately odd]], /(-V) = -/(V)- The 
sharp changes in the kick function correspond to close encounters with the planet, which occur at 




± 


arccos( 2 ? - 1 ) - - ?)'''^( 2 ? + l)/3 


(13) 


Kick functions for different perihelion distances are shown in Figure 1(b). (Everhart (1968) made an 
impressive early effort to obtain an empirical £t to the distribution of energy changes caused by planetary 
perturbations; however, his results are averaged over -0 and thus are not useful for our purposes.) 

In our mapping (|^ we represent /(V’) by a continuous interpolation function 


FW 


-F(-^) = 


IV' -V'ol E + 

n=l 

Clip + Co 

iv-- V’or'-'" E c<+V” 

n=l 


— TT < Ip < Ip-, 


lp-<lp < lp+, 


ip+<ip<0, 


(14) 


where we have divided the interval [—vr, 0 ] in three parts [—7i,'ip-], {ip-,ipjP), [' 0 +,O] and dehne ip-, - 0 +, 
and 00 by 

/(0_) = min /(0), /(0+) = inax /(0), 0o = (0_ + 0+)/2. (15) 

i^G[—7T,0\ 7r,0] 

Of course the function F(0) has the same odd symmetry as /(0). Note that only 10 of the 12 coefficients 
are independent, because of the continuity constraint at 0_ and 0+. 


The coefficients of the interpolation function F{'ip) are given in Table 1, along with the diffusion 
coefficient Dq dehned by equation (^) for initially parabolic orbits with the same inclination distribution. 


3. Results 

We examine a map with planet mass rrip = 5.24 x 10“®, which corresponds to the mass of Neptune, 
using the kick function for perihelion distance q = 0.5 (solid line in Figure 1(b)). 

Figure 2(a) illustrates the behavior of the map. The vertical axis is E/27i^ (the normalization is 
chosen so that the energy of an orbit with the planet’s semimajor axis is unity), the horizontal axis is 
0 / 7 r, and each point corresponds to one iteration of the map (one perihelion passage). The initial energy 
and azimuth of the orbits plotted are distributed uniformly random in the range —2 < E/2n^ < 0 and 
—TT < 0 < TT. The Figure shows the evolution of 100 comets over a time interval 300 (recall that the 
period of the planet orbit is unity). Figure 2(c) shows a magnified view of a smaller energy range, 
-1.50 < E/271^ < -1.28. 

Figures 2(b) and 2(d) are the same as Figures 2(a) and 2(c) except that the evolution of the comet 
orbits is determined by direct integration of the restricted three-body problem. The similarity of the 


^Strictly, the energy change is only odd to first order in rup. 
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Fig. 1.— (a) The points show the energy change AE{ip)/mp for a sample of initially parabolic comets 
with perihelion distance q = 0.5. The solid line is the kick function {AE{'ip))/rrip = flip) for this 
perihelion distance, (b) Kick functions for perihelion distances g = 0.1 (dashed line), q = 0.5 (solid 
line), and q = 0.9 (dot line). It is an artifact of the plotting program that the curves are not precisely 
odd in Ip near the spikes. 












Q 

0.1 

0.3 

0.5 

0.7 

0.9 

ijj-lTl 

-6.26 X 10-1 

-4.33 X 10-1 

-2.89 X 10-1 

-1.73 X 10-1 

-7.33 X 10-2 

xIj+Itt 

-6.23 X 10-1 

-4.29 X 10-1 

-2.86 X 10-1 

-1.69 X 10-1 

-7.00 X 10-2 

c[~^ 

6.56 X lOi 

1.47 X 10^ 

1.82 X 10^ 

1.77 X 10^ 

1.50 X 102 


7.26 

-4.08 

-6.22 X lOi 

-6.83 X lOi 

-5.37 X lOi 

C'f’ 

-1.76 X 10^ 

-9.64 X lOi 

-2.42 X lOi 

-9.93 

-1.07 X lOi 


5.56 X lOi 

1.76 X lOi 

2.01 X 10-1 

-1.15 

-2.35 X 10-1 


-2.40 X 10^ 

-2.14 X lO^ 

-3.41 X 10^ 

-6.48 X lO^ 

-2.45 X 10® 


-6.84 X lOi 

8.83 X lOi 

-6.60 X 10^ 

-2.94 X 10® 

-2.83 X lOi 


3.06 X lOi 

3.12 X lO^ 

-1.26 X 10^ 

-1.08 X lOi 

-2.39 X 10® 


1.04 X lOi 

1.62 X lO^ 

-8.11 X 10^ 

-1.22 X lOi 

-6.52 X 10® 


5.83 X 10® 

6.05 X 10^ 

3.91 X 10® 

2.82 X 10® 

1.44 X 10® 

C'l 

2.96 X 10® 

4.45 X 10^ 

4.32 X 10® 

5.24 X 10® 

6.46 X 10® 

Do 

3.9 X 102 

5.0 X 10^ 

7.7 X 10^ 

1.0 X 10® 

1.7 X 10® 


Table 1: The coefficients of the interpolation fnnction F{'ijj) and diffnsion coefficients (eq. §) for different 
perihelion distances. 
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corresponding plots suggests that the map captures most of the relevant dynamics in the restricted 
three-body problem. 

The hgures exhibit two main classes of orbit: (i) isolated islands consisting of regular orbits; (ii) a 
single connected stochastic orbit (the “stochastic sea”). The regular orbits are in resonance with the 
planetary orbital period. The future evolution of any point in the stochastic sea always terminates 
with escape {E > 0) or collision with the Sun (L < 0). The lower limit of the stochastic sea is at 
E = J = Einit — 2nLinit. As the energy approaches escape the fraction of phase space occupied by 
the stochastic sea grows. Notice that there are no KAM surfaces that extend across all values of the 
azimuthal angle V’; the reason is that orbits that suffer close encounters with the planet (eq. 0) are 
always stochastic. Thus stochastic orbits can escape from any energy range. Many of these features are 
reminiscent of the Fermi map with a non-differentiable forcing function (Lichtenberg and Lieberman 
1992). 

To explore the evolution of planet-crossing comets over much longer times we have used the map 
to follow 150, 000 comets with initial energy E/2tt‘^ = —0.242 (corresponding to semimajor axis 4.13 in 
units of the planetary semimajor axis) and initial azimuth 'ip distributed uniformly random in [0,27r). 
The initial semimajor axis was chosen so that all of the comets were in the stochastic sea and thus can 
escape. The comets are lost if they reach i? = 0 (escape) or Ej^E^ = —2.181 (collision with the Sun, 
as determined from (0)). We followed the comets for 2.7 x 10^ time units, which corresponds to the 
age of the solar system, 4.5 Gyr, if the planet has Neptune’s orbital period. After this time 3.7% of the 
comets survived. Of the comets that are lost, 96.9% escape and the remainder collide with the Sun. 
To show the phase portrait of the survivors, we plot every tenth perihelion passage for a time interval 
of 1000 Neptune years (Figure 3a). To magnify the detail, we also plot all of the perihelion passages in 
the interval —1.50 < E/2tt‘^ < —1.28 (Figure 3b). 

It is worthwhile to compare Figures 2(c) and 3(b). The most striking difference is that the islands 
of regular orbits are empty in 3(b); this of course is because the comets that are shown in 3(b) had 
initial conditions in the stochastic sea, whereas those in 2(c) were chosen randomly. A more interesting 
difference is that many of the points in Figure 3(b) are concentrated near the shores of the resonant 
islands (this is especially noticeable near E/2t:‘^ ~ —1.30, —1.41, —1.43 and —1.49). This is the well- 
known phenomenon of resonance sticking: stochastic orbits can be trapped for extended periods in the 
forest of tiny resonant islands near the shore of the stochastic sea (e.g. Karney 1983, Meiss 1992). We 
explore this phenomenon further at the end of this section. 


To examine the energy distribution of the surviving comets we divided the sample of 150,000 comets 
into two equal parts. In Figure 4(a) and 4(b) we plot the energy distribution of the surviving comets 
from each sub-sample after 4.5 Gyr; in particular we plot the normalized distribution of the number of 
perihelion passages per unit energy interval per unit time, i.e. the function 


/i(F;,f) 


n{E,t)/P{E) 

oo 

Jn{E,t)/P{E) dE 
0 


( 16 ) 


(see dehnitions in O- We see that the curves in the two panels have different shapes, which shows 
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Fig. 2.— (a) and (c) The perihelion passages of 100 comets over 300 planetary orbital periods (“Neptnne 
years”), as determined from the map ( ]T^) with perihelion distance q = 0.5. The comets are initially 
nniformly distribnted in energy and azimnth and the planet mass is rrip = 5.24 x 10“®. (b) and (d) The 
same as hgnres (a) and (c), compnted by direct nnmerical integration of the comet orbits. 
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Fig. 3.— The final phase portrait of 150, 000 comets, after 4.5 Gyr. The planet has Neptune’s mass 
and orbital period, and the comets initially lie in the stochastic sea, at El‘l'n‘^ = —0.242 (semimajor 
axis 4.13 in Neptune units), (a) Every tenth perihelion passage of all surviving comets, over an interval 
of 1000 Neptune years. The mapping structure for E/27r^ ^ — 1.8 is unreliable since very few comets 
penetrate to this region (for example, only 3% of the comets are lost through collision with the Sun at 
E= —2.181; all the rest are lost by escape io E > 0). (b) A magnified view showing all perihelion 
passages in the interval —1.50 < E/2t[‘^ < —1.28. 
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that most or all of the hne structure is statistical noise. 

Figures 5(a) and (b) show the depletion function N{t) as a function of time for several different 
initial energies. Each curve is based on an initial sample of 20, 000 comets. The survival fraction at 
large times in these plots is unrealistic, as it consists almost entirely of comets that have been kicked 
into near-escape orbits of very long period (in practice these comets would be removed by tidal forces). 
The expected asymptotic behavior as f > cxd is easy to derive heuristically. Let us imagine giving a set 
of bound comets a broad distribution of energy kicks at f = 0. Since the distribution of kicks is smooth, 
the number density of comets as a function of energy is flat near E = 0, n(E)dE = const x dE. The 
distribution of orbital periods is then n{P)dP = n{E)\dE/dP\dP = const x P~^^^dP. These comets 
remain bound until their second perihelion passage, after which they will normally be ejected within a 
relatively short time. Thus the expected asymptotic shape is 

OO 

N{t) oc J n{P) dP oc (17) 

t 

This power law is shown by an arrow on Figure 5(b) and accurately reproduces the behavior of the 
map. 

A more realistic map can be obtained by choosing a slightly negative escape energy, to crudely 
represent the effects of Galactic tidal forces (e.g. Heisler and Tremaine 1986). In Figures 6 (a) and 
(b) we plot N{t) as determined by mapping 500,000 comets with initial energy £'/27r^ = —0.242 and 
random initial phases. The two solid lines correspond to two different choices for the escape energy, 
£'esc/27r^ = 0 and —0.03; the latter corresponds to the more realistic case where a comet escapes from 
the influence of the giant planets when it reaches a semimajor axis a ~ lOOOAU where galactic tides 
and passing stars begin to have signihcant effects (Duncan, Quinn and Tremaine 1987). When the 
escape energy is non-zero N{t) decays approximately exponentially [Figure 6 (a)] until G ~ 5 x 10^ and 
thereafter as t~^, k ^ 1.3 [Figure 6 [b]). We interpret this late-time behaviour as the result of resonance 
sticking. 

We may also compare the depletion function in the map to the predictions of the diffusion approx¬ 
imation (eq. H). The dotted line shows the analytic prediction from equation (H), using the diffusion 
coefficient Dq from Table 1. As we indicated in § | 1 . 1 | , a more accurate approximation is that comets 
are lost either when E > 0 (escape) or E/2tt‘^ < J= —2.181 (collision with the Sun). We have 
also solved equation (H) numerically with these boundary conditions to provide a second prediction of 
N{t) (dashed hne). For reference, we show by solid arrows on Figure 6 (b) a power-law depletion law 
N{t) oc with k = 1.3, fc = | as predicted by equation (P!7|), and k = 2 as predicted by the diffusion 
approximation with Emin —^ —oo (eq. H). 

We close this section by exploring the phenomenon of resonance sticking in more detail. To do 
so, we map 5 x 10® comets with initial energy Ej2'n'^ = —0.49, and consider that a comet is lost 
anytime it reaches the nearby boundaries Eminl‘^'^‘^ = —0.59, Emax/2vr^ = —0.39 (focusing on a smaller 
energy interval provides a more sensitive probe of sticking, since the diffusion time is reduced for hxed 
diffusion coefficient). The resulting depletion functions for the mapping (eq. 0) and the diffusion 
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Fig. 4.— The energy distribution of the surviving comets from Figure 3, divided into two equal sub¬ 
samples. The energy distribution has been averaged over 10^ Neptune years to smooth out short-term 
fluctuations. The spikes are generally not common to the two plots and hence represent statistical noise. 
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Fig. 5.— The logarithm of the fraction of surviving comets in the map ([T^) as a function of (a) 
time and (b) log of the time. Each sample initially contained 20, 000 comets uniformly distributed in 
azimuth, with initial energies E/2t:‘^ = —0.985, E/27r^ = —0.490, E= —0.242 and £'/27r^ = —0.03. 
The initial energies are chosen so that all comets lie in the stochastic sea. The slope of the asymptotic 
law 0) is shown by an arrow. 
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Fig. 6.— The solid lines show the logarithm of the fraction of surviving comets in the map ( p!^ 
as a function of time (a) and logarithm of time (b). The sample initially contained 500,000 comets 
with energy i?/27r^ = —0.242; these were followed for two values of the escape energy, ii^esc/2vr^ = 0 
and F^esc/ 27 r^ = —0.03. We also plotted two predictions from the diffusion approximation of § |1 . 1| . The 
dotted line is for a single absorbing boundary at = 0 (escape) and the dashed line is for two absorbing 
boundaries at i? = 0 (escape) and = J/27r^ = —2.181. The asymptotic falloffs N{t) oc 

N{t) oc and N{t) oc t~^'^ are shown by arrows. 
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approximation (eq. ^ are shown in Figure 7 by solid and dashed lines respectively. We see that for 
t ^ 7x 10^ the depletion for the mapping is exponential, just as predicted by the diffusion approximation, 
although the rate of depletion in the map is slower than the diffusive prediction. However, for f ^ ~ 

10® the map again shows a power-law depletion, N{t) oc instead of the exponential depletion 

predicted by the diffusion approximation. Figure 8(a) shows the map over 1000 Neptune years for the 
67 comets that survived for t = 5 x 10®. Most of these comets are at the boundary of resonant islands, 
conhrming that the power-law tail of the depletion function is due to resonance sticking. 

Indirect evidence for resonance sticking is also seen in the energy correlation function p{AE,At), 
which gives the probability that a comet suffers an energy change AE over a time interval At. For 
\AE/E\ -C 1 the diffusion approximation predicts 


1 

2 _ DAt 

where P{E) is the comet’s orbital period and D is the mean-square energy change per perihelion 
passage (eq. ^). We have plotted equation (^) as a dotted line on Figure 8(b). We have also plotted 
the correlation function for the map, at time intervals At = 20, 60 and 100 (to estimate the correlation 
function we followed 10® comets with initial energy E/2tt‘^ = —0.49 for 10® Neptune years; we calculated 
the normalized distribution of energy changes after every At time units and averaged these distributions 
to obtain the correlation function). We see that as At increases, the correlation function for the map 
becomes more and more sharply peaked, which is a sign of resonance sticking. A possible concern is 
that our results are biased by a selection effect: we only followed comets so long as they remained 
inside a narrow energy band, —0.59 < i7/27r^ < —0.39. To check whether this bias is present we have 
constructed correlation functions using a Gaussian random walk with mean-square energy step given by 
the diffusion coefficient D and the same boundary conditions; these are shown in Figure 8(b) as dashed 
curves. In fact the dashed curves are so similar that they appear as a single near-Gaussian curve, which 
confirms that the peak seen in the correlation function for the map is not a result of the absorbing 
boundaries. 




exp 




(18) 


3.1. Other planet masses 

So far we have explored the Keplerian map for planet mass rup = 5.24 x 10“®, corresponding to 
Neptune. Figures 9(a) and (b) show depletion functions for Saturn (mp = 2.86 x 10“'^) and Jupiter 
{rrip = 9.55 x 10“^). In both cases we calculated time in the corresponding planet years and comets 
were lost when they either reached semimajor axis 1000 AU or collided with the Sun [L < 0). The 
Saturn and Jupiter maps began with 2 x 10® and 10^ comets respectively. The results for all three 
planet masses are similar: the number of surviving comets decays approximately exponentially until a 
time ts, and thereafter decays as a power law, N{t) oc with k ~ 1.3. The time tg appears to be 
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Time Logjg(Time) 

Fig. 7.— The solid lines show the logarithm of the fraction of surviving comets in the map (13) inside 
the energy interval —0.59 < £'/27r^ < —0.39, as a function of time (a) and logarithm of time (b). 
The results are based on an initial sample of 5 x 10® comets with initial energy E/2tt‘^ = —0.49. The 
depletion function for the diffusion approximation (^) is shown by a dashed line. The arrow shows the 
power law N{t) cx 




E/27T 
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Fig. 8.— (a) The map of the 67 comets that survived after 5 x 10® Neptune years of the mapping 
shown in Figure 8. The orbits are stuck to the edges of the resonant islands, (b) Correlation functions 
ap{AE/a) for time intervals At = 20, 60, 200 (solid lines). The Gaussian correlation function expected 
from the diffusion approximation (eq. and the correlation function for a Gaussian random walk 
with the same boundary conditions as the map are shown by dotted and dashed lines respectively. 
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independent of the initial energy: ~ 4 x 10^ planet years for Neptune (from Figure 6[b]), 1.7 x 10® 

for Saturn, and 1.5 x 10® for Jupiter. These times can be crudely £t by the formula tg ~ 0.14/mp. 

Figure 10(a) shows 1000 Jupiter years of the map, for the 47 comets with initial energy E/2tt'^ = 
—0.5 that survived for 2 x 10® Jupiter years. Once again, most of the surviving comets are stuck near 
resonant islands. 

To compare the depletion functions for Jupiter, Saturn and Neptune we superimpose these on 
Figure 10(b) for (approximately) the same initial energy ~ —0.24. Note that the timescale is 

planet years. 


4. Discussion 


The results of the previous section suggest that the depletion function for the Keplerian mapping, 
with absorbing boundaries at 0 < E^ax < E < E^ini can be approximately represented in the form 


N{t) ~ exp(—ampDof/P), t < tg 

~ fiitjtf, t > t,. 


( 20 ) 

( 21 ) 


Here the time t is measured in physical units (e.g. Earth years), P is the planet’s orbital period in 
Earth years, a and /3 are constants that are approximately independent of planet mass but may depend 
on the initial energy or angular momentum of the comets, Dq is dehned by equation (^), and k ~ 1.3. 
This formula only applies if the escape energy is non-zero. The “sticking time” tg is given by 


log(l//3) P 

a Do 


( 22 ) 


We show sticking times calculated by equation (|2^) , with a 
as vertical dotted lines. 


0.04 and /3 10 on Figures 9 and 10(b) 


Our studies of the Keplerian map show that resonance sticking does occur for highly eccentric 
planet-crossing orbits. The characteristic sticking time, after which resonance sticking dominates the 
depletion function, is — 2 x 10® (Earth) years for Jupiter, 5 x 10^ years for Saturn, 5 x 10® years for 
Uranus, and 7 x 10® years for Neptune (for perihelion distance 0.5 times the planet’s semimajor axis). 
Thus resonance sticking would be unimportant for a planetary system containing one planet with the 
mass and semimajor axis of Uranus or Neptune. However, if the planet resembled Jupiter or Saturn, 
resonance sticking would dominate the long-term behaviour of comets on planet-crossing orbits. 

These results invite comparison with long integrations of the orbits of Neptune-crossing objects. 
Duncan and Levison (1997) followed 2200 test particles on Neptune-crossing orbits for 4 x 10® years 
in the gravitational held of the Sun and the giant planets. They found that on Gyr timescales the 
depletion function could be described by a law of the form (^) with fc ~ 1, normalized so that 1% of 
their particles survive after 4 Gyr. 
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Logjo(Time) 



Fig. 9.— The depletion functions for the Keplerian map with a planet mass equal to that of Saturn 
(a) and Jupiter (b). We used different initial energies for the comets, as indicated in the Figures, but 
in all cases the comets were in the stochastic sea. The comets are assumed to escape when they reach 
a semimajor axis 1000 AU or collide with the Sun (T < 0). The vertical dotted lines represent sticking 
times dehned by equation 







E/27T' 
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Logjo(Time) 


Fig. 10.— (a) The map of the 47 comets that survived after 2 x 10® planet years of the mapping, 
with Jupiter’s mass and initial comet energy = —0.5. The orbits are stuck to the edges of the 

resonant islands, (b) A superposition of the depletion functions for mappings with Jupiter, Saturn and 
Neptune [from Figures 6(b), 9(a), 9(b)]. The initial comet energy (in the corresponding planet energy 
units) was approximately the same in each case, ~ —0.24. The vertical dotted lines represent 

sticking times dehned by equation (^^. 






























22 


These results are reminiscent of the results from the map: there is power-law behaviour at large 
times, the exponent A: ~ 1 is close to the exponent 1.3 observed in the map, and the surviving comets 
at the end of the integration were mostly in Neptune resonances, just as in the map. However, our map 
predicts that sticking is only important for t > tg = 7 x 10® years, and that aX t = tg only 10“^ of the 
initial particles survive, which is smaller by two orders of magnitude. The reason for this discrepancy 
between the map and direct integrations is unknown; perhaps the enhanced role of resonances reflects 
the extra degrees of freedom in the integration, which has four planets on eccentric, inclined orbits 
rather than one on a circular orbit. 

To investigate evolution in more realistic planetary systems, we would like to generalize the map 
to include (i) multiple planets; (ii) eccentric planetary orbits; (iii) dependence of the perturbing Hamil¬ 
tonian K on energy and angular momentum, not just the azimuthal angle 'ijj (cf. eq. 0). All of these 
generalizations are straightforward in principle, although the kick function will depend on more vari¬ 
ables and hence requires a more complicated parametrization. When k depends on E and L, integrating 
Hamilton’s equations across the delta-function at £ = 0 is generally no longer analytic. One possibility 
is to use a heuristic form for k chosen so that the integration across £ = 0 is analytic. A more general 
approach is to work with the mixed-variable generating function that describes the canonical transfor¬ 
mation induced by the perturbing Hamiltonian across i = 0. Thus if (L, —E, zu, t) are the canonical 
variables just before perihelion passage at £ = 0 , primes denote the same variables after perihelion 
passage, and S{L', E',zu,t) is the generating function, then 

. _ dS dS , _ dS , _dS 

It is also straightforward to show that 

S'(L', E', w, t) = L'vj — E't — mpK^L', E\ tu, i) + 0(mp). (24) 

Equations ( [23|) are implicit in the post-perihelion variables, but iterative schemes converge rapidly 
for small planet mass rup. 


(23) 


5. Summary 

The Keplerian map offers considerable insight into planet-induced evolution of highly eccentric 
orbits, both because of its simplicity and its computational speed. 

The map reveals that the exponential decay in the number of planet-crossing comets predicted by 
the diffusion or random-walk approximation is replaced at late times by a power-law decay, N ~ 
and that most surviving comets are “stuck” at the edges of resonant islands. Resonance sticking 
dominates the evolution of comets under Jupiter’s influence after timescales as short as a few Myr. 

Sticking to Neptune resonances is not strong enough to explain the remarkably high survival fraction 
of Neptune-crossing comets (1% after 4 Gyr) found by Duncan and Levison (1997), most likely because 
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the simple map used in this paper does not include eccentric planetary orbits, multiple planets, or 
dependence of the planetary perturbations on perihelion distance. 

An unresolved issue is how small perturbations and dissipative effects (passing stars, non-gravita- 
tional forces, etc.) affect the lifetimes of comets that are stuck to resonances. 

We thank Martin Duncan, Matt Holman and Norm Murray for discussions. This research was 
supported in part by NASA grant NAG5-7310. 
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